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We propose a three-component reaction-diffusion system yielding an asymptotic 
logarithmic time-dependence for a moving interface. This is naturally related to 
a Stefan-problem for which both one-sided Dirichlet-type and von Neumann-type 
boundary conditions are considered. We integrate the dependence of the interface 
motion on diffusion and reaction parameters and we observe a change from transport 
behavior and interface motion ~ t^/^ to logarithmic behavior ~ In t as a function of 
time. We apply our theoretical findings to the propagation of carbon depletion in 
porous dielectrics exposed to a low temperature plasma. This diffusion saturation is 
reached after about 1 minute in typical experimental situations of plasma damage 
in microelectronic fabrication. We predict the general dependencies on porosity and 
reaction rates. 

I. INTRODUCTION 

The generation and propagation of interfaces in react ion- diffusion systems is of great 
importance in a variety of scientific and technological contexts. In nature, the interface often 
spatially separates two macroscopically different homogeneous regimes in which one phase 
or product dominates another one. Typically, at least one particle-type diffuses through 
the medium while reacting with other particles. Births and deaths depend on the spatial 
location and on the presence of other particle-types. Then, when installing appropriate 
initial and boundary conditions one often sees the appearance of an interface and one asks 
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for the typical temporal behavior in its initial growth and further propagation. The case we 
treat here is boundary driven from one side, with symmetric bulk diffusion for one particle 
type that invades the region of a second particle type, while constantly being annihilated by 
yet a third type of particles. That creates a moving phase boundary and thus goes under 
the general name of Stefan problem jl, 2]. 

The present work is inspired by experimental investigations on plasma damage in porous 
low-capacitive materials (low-Zc, see j3|). One imagines here a porous SiOC material which 
is exposed via one of its physisorbedaries to radicals such as active oxygen atoms. These 
radicals diffuse through the material while removing hydrophobic carbon-hydrogen groups 
and while being chemi-adsorbed in the pores of the material. One then observes the appear- 
ance and motion of an interface separating a rich (pristine) carbon-phase from a damaged 
(oxygen-rich) layer. It is important here to understand the time-dependence of that motion 
and its dependence on diffusion and reaction parameters, which can further be related to 
material and chemical properties. We come back to this application in the separate section 

Of course, the literature on interface motion from reaction-diffusion equations is vast, and 
many different phenomena have been reported. In particular, damage spreading and related 
front propagation can take many different forms which we cannot at all review here; see 
e.g. for general research aspects. As just one very famous example for one-dimensional 
front propagation, as we also have here, react ion- diffusion equations such as of the type of 
Kolmogorov-Petrovsky-Piskounov show interface motion between two homogeneous phases 
at constant speed That then takes the form of travelling waves in which one phase 
advances over the other phase. Here, however, the proposed model shows the motion of the 
interface x{t) ~ ln(l + at + b\/t) (in one dimension as function of time t) to be logarithmic 
which means that after an initial period diffusion saturation occurs. The main point is that 
the diffusive motion is boundary driven and that the propagation is limited by the reactions. 
For the purpose of the application, the logarithmic behavior is caused by the disappearance 
of the oxygen radicals before reaching the moving interface separating them for the carbon- 
rich environment. The -y/t within the logarithm, mostly visible in an initial time-period, 
originates directly from the diffusive motion. 

The next section contains the specific mathematical model in the form of coupled partial 
differential equations. They give a mean-field description of the react ion- diffusion system. 
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The main finding is the change from a diffusive to a specific logarithmic time-dependence 
for the interface x{t) as above, with the study of dependencies on initial conditions 
and dynamical parameters. We give a physical heuristics for understanding the basic 
characteristics via a linearized model. 



The section Uni on numerical results summarizes the behavior of our model, as a function 
of the physical parameters. Our results agree with the analytical results of the linearized 
model, but the higher flexibility of simulations allow for obtaining additional information 
and predictions concerning the dependence on the various parameters in the model. 

Section [IV] goes deeper into the specific application that triggered our study, that of 
the propagation of damage in porous low-k materials. We recall the main aspects and the 
meaning of our result for that context. In particular, reference js] contains already the basic 
set of differential equations, but no systematic analysis was performed and the distinction 
between different (one-sided) boundary conditions was not considered. The present paper 
also adds the interesting connections with reaction-diffusion systems of Stefan type. 



II. MODEL AND MAIN FINDINGS 



There are three types of particles. The first are called free 0-particles, which can diffuse 
but can also change into bound 0-particles and further annihilate with either bound O- 
particles or with C-particles: 

free O 4- bound O, free O + bound O A 0, free O + C 4- 

and no other reactions or diffusions take place. The letters are chosen in the context of the 
later application in Section IIVI with O for oxygen and C for carbon. The word "free" will 
there be understood as "reactive" and "bound" will mean adsorbed to the inner surface, 
either physically or chemically. Once oxygen atoms are absorbed on the surface they start 
to move around until they recombine and desorb. In other words, since in the experimental 
conditions desorption of bounded oxygen molecules typically occurs through recombination, 
reversible desorption is not taken into account, jo], |lo|. For the moment we continue to use 
the more general words "free" and "bound," as the mathematical set-up can have different 
natural realizations. 
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The dynamics is further determined by the initial condition at time zero when there are 
only C-particles with homogeneous density, and by the boundary conditions. The physics 
we have in mind concerns a three-dimensional region of locations (x > 0,y,z), which is 
isotropic over the (y, 2;)— coordinates. For the boundary conditions at the x = plane we 
consider two possibilities. The first, of Dirichlet-type, specifies a constant density of free 
0-particles at x = 0; the second, of von Neumann-type, gives a constant flux of particles 
entering the region at x = 0. For the questions and the problem at hand we restrict ourselves 
to an effective mean field treatment, ignoring fluctuations or microscopic inhomogeneities, 
and the geometry is one-dimensional. 

We consider the coordinate x > in which to express the profile of particle densities as a 
function of time t > 0. There are thus three types of concentrations p(x, t), c(x, t) and m(x, t) 
respectively for free 0-particles, C-particles and bound 0-particles. The reaction-diffusion 
system is then defined via the equations 

p = D p" — K pm — Rc p — a p 
c = —Rep 

m = —upm + ap (H-l) 

where the dotted left-hand sides are time-derivatives. The concentration p (of the free O- 
particles) undergoes a second spatial derivative in the right-hand side of the first line of 
flll.ip : the parameter D > being the diffusion constant. The other parameters k, R and a 
are positive reaction rates. The initial condition {t = 0) specifies 

c(x, t = 0) = 1, p(x, t = 0) = m(x, t = 0) = 0, for X > (II.2) 

meaning that we start with C-particles fully occupying the material. On the other hand, 
the boundary conditions (x = 0) can be taken either 

Dirichlet-type: p(x = 0,t) = V, or von Neumann-type: p'(x = 0,t) = — J (II-3) 

where we can take V = 1 without loss of generality; the J > stands for the rate of free 
0-particles entering the region which for the moment we leave as a parameter. It seems 
more natural to opt for von Neumann-type boundary conditions for the specific application 
we have in mind, but we checked that after a transient time there is really no difference 
between Dirichlet and von Neumann conditions. In the case of von Neumann-type boundary 
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conditions, the value p{x = 0, t) converges from being zero at time zero to a fixed finite 
value (equal to some V > depending on J) exponentially fast as time evolves. As a 
consequence, we mostly concentrate on just one type, that of Dirichlet boundary conditions 
with p{x = 0, t) = 1 for all times, but further discussion on possible differences will follow 
below. 

The results of numerical integration of the differential equations (lll.ip are detailed in 
Section IIIII Our main finding is that the concentration c{x, t) of C-particles (defining the 
absence of damage in the application of Section irv]) shows an interface x{t) separating a 
C-poor from a C-rich phase. There are a priori a number of different possible mathematical 
definitions of the position of the interface x{t). A first one is to solve c{x{t),t) ^1/2 
(or some other number) for x(t), which certainly makes sense when the interface is sharp 
enough, after some transient time, but is less useful for very small times. A second definition 
uses the same conditions as above but we add the equation x{t) = —Dp'{x{t),t) for the 
interface speed. That brings the problem in the natural neighborhood of Stefan problems 
ll, |2| , also valid for very small times but more involved numerically. Whatever (reasonable) 
definition we take for x{t) we invariably find that the motion is logarithmic of the general 
form x{t) = 741n(l + at + by/t) in a typical range of parameters for diffusion and reaction 
rates k, R and a. The coefficients A, a and h depend on these parameters, see also in 
Section [Inl For small times the motion is diffusive x{t) ~ \/t while it saturates as x{t) ~ Int 
for large times t. The transition between the two regimes is basically determined by the 
amount of reactivity, the more the 0-particles can be annihilated, the faster saturation sets 
in. 

An approximate treatment towards the second Stefan-like definition of the interface mo- 
tion described above can already be made visible for a much simpler linear system to which 
we turn next. We simplify (III.ip to a one-component react ion- diffusion equation by imagin- 
ing an autonomous linear dynamics for p of the form 

p{x,t) = Dp" -Pp (II.4) 



The equation ( III.4P mimics the first line of in the approximation nm + Rc+a = /3 > 



constant, which we expect to be mostly reasonable when c = and m =constant which 
makes sense in the damaged region (0 < x < x{t)) where the carbon was depleted. The 
right boundary denoted by £ in what follows therefore plays the role of interface. We also 
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can think of /3 as a measure of reactivity and interaction, while D measures the diffusion. 
The main advantage is of course that now the equation flll.4p can be solved exactly with 
boundary conditions p{x = 0,t) = l,p{x = i,t) =0 and initially p{x,t = 0) = 0. The 
solution is 

p[x,t) = pi[x) + e ^ 2_^(^nG sm — X, where 

1 

nvr 



l-e2fe« 1-e 

oo 



= — ^^a„sin— X which determines the a„ (H-S) 

1 ^ 

and p^ is the stationary solution on [0,£]. In fact we see that p{x,t) converges exponentially 
fast to p£ with rate /3 > uniformly in i. We can already learn a great deal from the 
behavior at the edge x = i. It is easy to find that 

We think now of i as the position x of the interface and —p'i{i) is then the speed of the 
interface, so that Stefan's condition x{t) = Dk/ sinh(A;x(t)) can be integrated to give 

x{t) = ln[l + /3t + ^/2^t^^] (II.7) 

We clearly recognize the behavior x{t) oc ^J~DJP> Int for large times t resulting in a straight 
line in Fig. 1 and intersecting the vertical axis (small t) at a/-D//3 ln(2/3). For small times t 
flll.Tp gives X (t) oc Vt, and x{t = 1) ~ V2D for small /3 > 0. We can thus find the effective 



D and /3 for (111.4^ from the intersections of the x(t) as function of Int with the vertical 
axis, as in Fig. 1, which provide consistency checks for our modeling equations. 

Similarly, we can also solve flll.4p with von Neumann- type boundary condition p'{x = 
0, t) = — J but still putting p{x = i,t) = and initially p{x, t = 0) = 0. The solution is 

/ \ / \ -Bf'^, _D_(n-i/2)f^^ (n — l/2)7r 
p[x,tj = ps[xj + e ^ 2,^ne ^ cos x, where 

1 ^ 

^ (ti 1 /SItt 
= — bn COS X which determines the bn- (H-S) 
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Exponentially fast in time t, p(0,t) = J/k + 0{e ^^), and hence we find that for J ~ 




the von Neumann condition imitates the Dirichlet condition with fixed p(0, t) = 1 and 
p^{x) ^ Ps(x) ^ exp —kx for large at least for (3^0. Clearly then, when the rate /3 would 
go to zero, some differences between Dirichlet and von Neumann (one-sided) boundary 
conditions can be expected, as we indeed recover in Section IIII Bl for the real model. The 
next section goes back to the full nonlinear model equations (lll.ip but we will find the many 
general features to be similar to our chosen approximation above. 



III. NUMERICAL RESULTS 

We now turn to the full set of equations (lll.ip . and switch to computer simulations. First, 
we study the system for a specific set of parameters D = 1.0, k, = 0.1, R = 1.0 and a = 0.1, 
and apply Dirichlet boundary conditions; the primary motivation for this specific choice of 
the parameters lies in their numerical convenience. We apply a forward Euler integration 
scheme to equations (lll.ip with a time step of At = 0.001, and store the density profiles of 
C-particles, free and bound 0-particles at times t = 1,2,.... Based on our analytic solution 
flll.7p of the simplified description fill. 41) . we expect that for long times the interface position 
x{t), defined as the (linearly interpolated) location where the C-particle density equals 1/2, 
changes in time as x{t) ^ Aln{l + at), for some values of A and a which are determined by 
the system parameters D,k,R and a. Fig. 1 shows that this expected long time behavior is 
confirmed by the simulation data. 

In the same simulations, we also determined the concentrations of C-particles and free 
0-particles as a function of position, for the times t/100 = 1, 2, 4, 8, . . . , 8192. These are 
c{x,t) respectively p{x,t) in equations flll.l|) . The resulting profiles of C-concentration are 



shown in Fig. 2a. The C-density profile attains a sigmoid-like shape, which gradually moves 
deeper into the material. As the times of the curves are separated by a factor of two, and 
the depth increases with the logarithm of time, the curves are more or less equidistant. 

The concentration profile of free 0-particles is shown in Fig. 2b. The general shape is that 
of an exponentially decaying curve up to the interface position x{t), after which it decays 
again exponentially, but with a much steeper slope. Once the interface has significantly 
moved into the material, the gradient of free 0-concentration at the surface becomes constant 
even though we applied Dirichlet boundary conditions. This confirms our earlier analytic 
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FIG. 1: Position of the interface x{t) as a function of time. For longer times, the curve approaches 
a straight hne, indicating that asymptoticahy x{t) = Aln{l + at). 

observation that the behavior for Dirichlet and von Neumann conditions soon becomes 
indistinguishable. 

A. Influence of diffusion 

An important ingredient is to understand the influence of the diffusion constant D on 
the temporal behavior of the interface. After all, the porosity and the architecture of the 
material can modify exactly that diffusive behavior. 

We repeated the simulations with different values for the diffusion coefficient: D = 0.5, 1, 2 
and 4. Fig. 3 shows the rescaled interface position x{t) / a/D, indicating that the main effect 
of variations in D is a rescaled amplitude: A ~ V^, as expected from the linear model 
in fill. 71) . It appears therefore that we understand the influence of the connectivity of the 
porous material on the time-dependence of the interface. 

B. Influence of reactivities 

The adsorption of the free 0-particles and the chemical annihilation of the free O- 
particles with the bound 0-particles (e.g. to create a O2 molecule which is chemically 
neutral for C-particles) and with the C-particles is the reactive content of our model. 
Obviously the C-particles must remain to define an interface at all, but we can imagine 
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FIG. 2: Profile of tlie C-concentration (top panel) and free 0-concentration (lower panel) as a 
function of depth, at times t/100 = 1, 2, 4, ... , 8192. Note that the top panel has a linear vertical 
scale, while the lower panel has a logarithmic vertical scale. 

lowering very much the adsorption of the free 0-particles. In the extreme case we exclude 
the adsorption of free 0-particles by setting a = 0. Then, we have m(x, t) = always 
and everywhere which simplifies f lll.l|) somewhat. In the case of Dirichlet type-boundary 
condition p{x = 0,t) = 1, we basically have transport of 0-particles through a layer of size 
x{t) with difference in concentration fixed (and equal to 1): we thus expect from Pick's law 
that the current x{t) ~ ^/x(t) is proportional to the concentration gradient which gives us 
x{t) ~ We have checked that this behavior is indeed found for very small adsorption 
rates, a = and k = 0. Fig. 4 shows the numerical evidence. 
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FIG. 3: Rescaled position of the interface x{t)/\fD as a function of time, for simulations in 
which D = 0.5,1,2 and 4. For longer times, the curves fall on top of each other, indicating 
that asymptotically A oc y/D. 

On the other hand, with von Neumann boundary condition p'{x = 0,t) = — J and for 
very small absorption rates a = and k = 0, there is a continuous influx of free 0-particles 
leading to a concentration of them at the surface which increases in a power law fashion. 
That can be understood from the following calculation. Let us take J = 1 and consider the 
function 

satisfying indeed dtpi{x = 0,t) = —1. Then, dtpi{x,t) = pi{x,t)/t [1 +x/\/At\ and 
dlpi{x, t) = pi{x, t)/t, so that pi appears as a solution of the diffusion equation for x <^ y/i. 
We thus conclude that pi(0, t) grows as \/t. Numerically, we find that the lack of saturation 
of free 0-concentration at the exposed surface does not alter the asymptotic behavior of 
the interface, but it does lead to a logarithmic correction. The behavior is consistent with 
x{t) ~ -y/tlnt. Also for von Neumann boundary conditions, numerical evidence is presented 
in Fig. 4. 

IV. APPLICATION 

Low dielectric constant (low-k) materials are typically used in semiconductor manufac- 
turing to compensate for the RC-delay and power consumption increase in microelectronic 



11 



1000 - 



100 r 



10 r 



1 ' ' 1 ' ' ' 1 ' 1 

_ Uiricniet 


1 1 1 II 


: von Neumann 




- sqrt (x) 




sqrt (x*log(x)) ■^,< 




y^^ ' ' 

<^'/ 

- 7^/ 

- ' ^ 1 lJ 1 1 lJ 1 1 lJ 1 1 





1 10 100 1000 10000 100000 

t 

FIG. 4: Position of the interface x{t) as a function of time, for simulations in which the adsorption 
of free oxygens is excluded. With Dirichlet boundary conditions, the asymptotic behavior can be 



well fitted by x{f) ~ with von Neumann boundary conditions, this becomes x (t) ~ y/t Int. 

devices associated with their continuing scaling-down js], The microelectronics commu- 
nity has adopted /c as a representation of the relative permittivity or dielectric constant in 
contrast to the scientific communities using e^. Nowadays, porous low-k dielectrics based on 
silica and silsesquioxanes with 10-15% of organic hydrophobic groups are the most favoured 



class of materials for advanced interconnect technology nodes [Uj . The hydrophobic groups 
are bonded to a Si- atom in the Si04/2 matrix and can be represented as O3 =SiCHa.. To 
reach a /c-value below 3.0, introduction of artificial porosity is needed (the /c-value of air 
is close to 1) js]. Advanced \ow-k materials have a porosity of 25-50% and a typical pore 
size being close to 2-2.5 nm. Although the matrix of these materials has properties similar 
to traditional Si02, their chemical stability and reactivity strongly depend on porosity. 
Materials with interconnected pores are chemically active because of the high diffusion rate 
of active chemicals. In particular the surface diffusion of adsorbed radicals is much faster 
than the desorption, which is our motivation not to include reversible adsorption- desorption 



in our model equations (III.ip . see 
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for specific experimental estimates. 



During device fabrication, the porous \ow-k dielectrics are exposed to various plasmas 
that might degrade their properties. The most challenging plasma treatment is related 
to plasmas based on oxygen- or hydrogen containing chemistries that are used to remove 
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organic photo- resist masks used for pattern transfer from the \ow-k films [l^. The challenge 
arises from the similar chemical composition of the low-A; film and the organic photo-resist 
mask. Both materials contain C-H bonds. The plasma radicals that remove the organic 
mask also penetrate the porous low-k dielectric film destroying its Si-CHs bonds. This 
leads to formation of polar Si-OH groups as a result of direct chemical reaction with active 
radicals or saturation of =Si« broken bonds by OH groups from ambient moisture. The 



Si-OH groups are centers for further moisture adsorption [15|. The absorbed moisture with 
a /c-value of 80 can fill the whole open pore volume of the material significantly increasing 
its /c- value which is summarized as damage. 

Two commonly known approaches are used for the photo-resist mask removal: (i) a 
low temperature, low pressure anisotropic plasma, where the photo-resist is removed by an 
ion-assisted process, ox,d,.„g o. .educing p.a..a che,n.t*. at ,ow temperature. Q and 
(ii) hydrogen-based downstream plasma (DSP) where the resist is removed at high tem- 
)eratures by a thermally activated chemical process jioj]. According to recent publications 
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16|, both approaches remove C from low-Zc materials. However, ion-assisted processes 
remove hydrophobic Si-CHs groups and graphitized-C (residual-C that remains after film 
fabrication), while a thermally activated chemical process with hydrogen plasma removes 
only graphitized-C from low-A; dielectrics. In both cases the depth of C-carbon removal 
can be directly measured with spectroscopic ellipsometry js], [mI. Various experimental 
results have shown the typical time-dependence of the depth of C-removal. In particular 
the typical logarithmic time-dependence was observed, as the one we have derived in the 
present paper. 

The depth of the plasma radical penetration in the low-k pores determines the extent of 
plasma damage. Clearly then, the role of Knudsen diffusion is important and our analysis 
in Section IIII Al clarifies some aspects. For example, the extent of plasma damage can be 
a few times higher for mesoporous dielectrics with high porosity levels in comparison with 



nanoporous dielectric with low porosity levels 
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As a word of caution and limitations of our approach, plasma damage of low-k materials 
is a complex phenomenon that results in changes of their bonding structure and pore mor- 
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phology jst- The complexity also increases due to the fact that porous \ow-k dielectrics are 
usually amorphous materials with a random pore structure j^. Our set of reaction-diffusion 
equations flll.l|) can be generalized to include some of these aspects but soon a more detailed 
simulation analysis becomes necessary. 

V. CONCLUSIONS 

The propagation of damage when a pristine material is exposed to diffusing and reacting 
particles saturates logarithmically. The initial position of the interface is diffusive ~ \/t 
but soon saturates ~ Int when reactions prohibit the propagation of the damage. The 
problem can be modeled as a reaction-diffusion system of coupled differential equations 
where the position of the interface can be determined by viewing it as a Stefan problem. 
We have discussed Dirichlet and von Neumann boundary conditions, solved a linear 
approximation and we have found the detailed influence of reactivities and diffusivity. 
Numerical integration of the coupled reaction-diffusion equations is in full accord with 
experimental findings in the context of microelectronic fabrication of low-/c materials. 
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